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Abstract 

We show that the emergence of different surface patterns (ripples, dots) can be well understood by a suitable 
mapping onto the simplest nonequilibrium lattice gases and cellular automata. Using this efficient approach 
difficult, unanswered questions of surface growth and its scaling can be studied. The mapping onto binary 
variables facilitates effective simulations and enables one to consider very large system sizes. We have 
confirmed that the fundamental Kardar-Parisi-Zhang (KPZ) universality class is stable against a competing 
roughening diffusion, while a strong smoothing diffusion leads to logarithmic growth, a mean-field type 
behavior in two dimensions. The model can also describe anisotropic surface diffusion processes effectively. 
By analyzing the time-dependent structure factor we give numerical estimates for the wavelength coarsening 
behavior. 



1. Introduction 

In nanotechnologies, patterns of large areas are 
to be assembled in a cost effective way, possibly 
by self-organization processes. Different materials 
under different conditions have been shown to ex- 
hibit dot and ripple formation in a rather univer- 
sal manner [lj. Universality appears, when micro- 
scopic details of interactions become irrelevant, i.e. 
the diverging correlation length and the complex 
behavior of the system dominates. This is typical 
in nonequilibrium systems with currents 0, yj . It 
was pointed out 0, [B| that grooved surfaces and 
growth instabilities may emerge as the consequence 
of broken detailed balance condition 

P(C)R c ^c ± P(C')R c ^c (1) 

where P(C) denotes the probability of the state C 
and Rc-+c is the transition rate between states C 
and C . This means that complex structures and 
patterns can emerge during relaxation in nonequi- 
librium system. In case of ion-beam-induced mod- 
ification of surfaces the first theory to explain pat- 
tern formation was suggested long time ago Q . 

To demonstrate the general tendency of self- 
organization we present here simple case studies of 
fundamental models of statistical physics, showing 
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universal scaling and pattern formation. Under- 
standing of basic surface growth models still lacks 
many details, although the basic universality classes 
and important models have been introduced @, [1] • 
Analytical tools like continuum field methods have 
limited applicability, while numerical simulations 
have achieved limited precision. 

One of the simplest nonlinear stochastic differen- 
tial equation set up by Kardar, Parisi and Zhang 
(KPZ) describes the dynamics of surface growth 
processes in the thermodynamic limit. It specifies 
the evolution of the height function h(x, t) in the 
d-dimensional space 

d t /i(x, t) = v + ctV 2 /i(x, t) + A(V/j(x, t)) 2 + r?(x, t) . 

(2) 

Here v and A are the amplitudes of the mean 
and local growth velocity, a can be understood 
as a coefficient of surface tension driven smooth- 
ing and r\ roughens the surface by a zero- 
average, Gaussian noise field exhibiting the vari- 
ance: (?7(x,i)r/(x',t')} = 2L><5 d (x-x')(t-t')- The 
notation D is used for the noise amplitude and () 
means the distribution average. The surface scal- 
ing exponents are known in (1 + 1) d [lfj| . but for 
the case considered here, 2-dimensional surfaces of 
3-dimensional system, as well as in higher dimen- 
sions approximations are available only. The exis- 
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tence of a finite upper critical dimension, where a 
smooth mean-field behavior would enter is also de- 
bated. As the result of the competition of roughen- 
ing and smoothing terms, models described by the 
KPZ equation exhibit a roughening phase transition 
between a weak-coupling regime (A < A c ), governed 
by the Edwards- Wilkinson (EW) fixed point at 
A = [ll[ , and a strong coupling phase. The strong 
coupling fixed point is inaccessible by the perturba- 
tive renormalization group method. Therefore, the 
KPZ phase space has been the subject of contro- 
versies for a long time. 

Mapping of surface growth onto reaction- 
diffusion system allows effective numerical simula- 
tions and understanding of basic universality classes 
12, [ljj. In one dimension a discrete, restricted 
solid on solid (RSOS) realization of the KPZ growth 
is equivalent to the Asymmetric Simple Exclusion 
Process of particles [1J], while some of us have 
shown that the roof-top model mapping [TBI, can 
be generalized to higher dimensions [171 . |18| . 

In two dimensions one can map KPZ processes 
onto anisotropic, but oriented migration of directed 
dimers (l7l | . This mapping is interesting not concep- 
tually only, linking nonequilibrium surface growth 
with the dynamics of driven lattice gases 0, EH , but 
provides an efficient tool for investigating debated 
and unresolved problems numerically [18] . The sur- 
face built up from the octahedra can be represented 
by the edges meeting in the up/down middle ver- 
texes (see Fig. [1]) . The up edges in the x = x or 
X = y directions at the lattice site i,j are repre- 
sented by a x {i,j) — +1, while the down ones by 
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octahedron site deposition flips four edges, meaning 
two '+1 — 1' exchanges: one in the x and one in 
the y direction. This can be described by the gen- 
eralized update rule: 



(3) 



with probability p for attachment and probability q 
for detachment. We can also call '+l'-s as particles 
and '— l'-s as holes of a reaction-diffusion model on 
a square lattice. Thus, an attachment/detachment 
update corresponds to a single step motion of an 
oriented dimer in the bisectrix direction of the x and 
y axes. We update the neighborhood of the sub- 
lattice points, which are the crossing-points of the 
dashed lines. In 
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18[ we derived how this map- 
ping connects the microscopic model to the KPZ 
equation with A — 2^^ — 1 and investigated its 




Figure 1: Mapping of the 2 + 1 dimensional surface growth 
onto the 2d particle model (bullets). Surface attachment 
(with probability p) and detachment (with probability q) 
corresponds to exchanges of particles, or to anisotropic dif- 
fusion of dimers in the bisectrix direction of the x and y axes. 
The crossing points of dashed lines show the base sub-lattice 
to be updated. Thick solid line on the surface shows the y 
cross-section as a reminder for the one-dimensional roof-top 
model. When the shown desorption/absorption steps are ex- 
ecuted simultaneously they realize a surface diffusion jump 
of distance s = 3 along the y axis. 



surface scaling numerically. In particular we have 
shown that for the A = case our simulation results 
agree with those of the exactly known EW class. 

Surface diffusion is a much studied basic pro- 
cess |20j . Several atomistic models have been con- 
structed and investigated with the aim of realizing 
Mullins-Herring (MH) diffusion [H [22| and scal- 
ing (for a recent review see [U 



The Langevin 
equation of MH diffusion is a linear one, with a V 4 
lowest order gradient term 



d t h{x, t) = -kV 4 /i(x, t) + r?(x, t) 



(4) 



resulting from a curvature driven surface current. 
This equation is exactly solvable and exhibits scale 
invariancc of the roughness. In Sect. 11.21 we show 
how this can be realized using the octahedron 
model. 

1.1. The simulations 

Dynamic, bit-coded simulations were run on con- 
served lattice gas models of sizes LxL following the 
rule ©. The surface heights were reconstructed 
from the slopes 



(5) 
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The interface width defined as 



w(L,t) = 



L 



1/2 



(6) 

has been calculated at certain sampling times t. 
In the absence of any characteristic length, growth 
processes are expected to follow a power-law behav- 
ior and the average surface width W — w can be 
described by the Family-Vicsek [13] scaling law: 

W(L,t)^L a f(t/L z ), (7) 

with the universal scaling function /(it) of the form: 

•P if u<l 
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const. 



if u > 1 
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Here a is the roughness exponent of the stationary 
regime, when the correlation length has exceeded 
the system size L and f3 is the growth exponent, 
describing the intermediate time behavior. The dy- 
namical exponent z is the ratio 

z = a/p. (9) 

For the 2d KPZ universality class our numerical 
results [H| 

a = 0.39(1), 



(3 = 0.245(5), z = 1.6(1) (10) 



are in agreement with most of the literature values 
HI, 2(|. The MH class, with non-conservative noise 
is characterized by the exponents 



a=l, 13 = 1/4, z = 4 



(11) 



in two dimensions, while in case of conservative (dif- 
fusive) noise the the exponents are 

a = 0, 13 = 0, z = A (12) 

and the growth can be logarithmic Q. Note that 
in both cases the the dynamics is very slow, de- 
scribed by the large dynamical exponent z = 4, 
which makes the simulations time consuming. 

Emerging patterns can be characterized by the 
time dependent power spectrum density (PSD) of 
the interface 



S(M) = (fc(k,t)M-M)> , 



(13) 



where the height in the Fourier space is computed 
as 



h(k,t) 



1 



L d/2 



£>(X,«)-S|exp(*kx). (14) 



We computed h(k,t) from the surface profiles with 
the FFT method and determined 5(k, t) by aver- 
aging over x and y directions in case of isotropic 
patterns and only over the direction perpendicular 
to the ripples in case of x/y anisotropy. We have 
also calculated the characteristic wavelength of the 
patterns from the maxima of S(k, t). 

The simulations were started from a flat surface, 
corresponding to a zig-zag configuration (alternat- 
ing T' and '0' lines of the lattice gas) of the slopes 
(see FigfT]) with periodic boundary conditions. In 
the model each lattice site can be characterized by 
16 different local slope configurations, but we up- 
date them only when condition ([3]) is satisfied. Fur- 
thermore, due to the surface continuity, two slopes 
become redundant, when considering the neighbors. 
Thus we can describe a lattice site by using only 
two bits. This permits efficient storage manage- 
ment and large system sizes L. The updates can 
be performed by bit operations either on multiple 
samples at once or on multiple sites in parallel. Our 
multi-sample algorithm proved to be ~ 40 times 
faster than the conventional FORTRAN 90 code, 
while a speedup factor of ~ 240 was measured us- 
ing the multi-site CUDA code on a NVIDIA C2070 
graphics card in comparison with a single 15 CPU 
core running on 2.8 GHz. The latter results could 
be achieved on huge system sizes up to L = 2 17 . 
More details are published elsewhere 17, 27, 28j |. 



1.2. Generalizations of the octahedron model 

An obvious first step of generalization is to com- 
bine the deposition and the removal processes, cre- 
ating a conserved dynamics. A simultaneous oc- 
tahedron detachment and deposition in the neigh- 
borhood can realize an elementary diffusion step 
(see Fig. Q). In ^ we described RSOS models 
with Ah = ±1 height restriction, which in the limit 
of weak external noise exhibit MH scaling. Micro- 
scopic models realizing this behavior are mainly un- 
restricted solid on solid (SOS) type, which can pro- 
vide steep slopes and strong curvatures. However, 
it turned out that the asymptotic universality class 
of the various limited mobility growth models is a 
surprisingly subtle issue and some earlier findings 
proved to be incorrect due to pathologically slow 
crossover and extremely long transient effects. 

In our model the target site is chosen in the ±x 
or ±?/ direction, with the probabilities: D +x , D- x 
or D +y , D- y respectively (see FigfT]), thus we can 
follow anisotropic surface diffusion in principle. 
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Throughout of our studies we normalized the at- 
tempt probabilities. The maximal jump distance 
was fixed to be s m < 4 lattice units following com- 
puter experiments. To control the surface diffusion 
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Figure 2: Local slope a x (i,j) (circles) and curvature c x (i,j) 
(squares) variables at an update site. Filled circles corre- 
spond to upward, empty ones to downward slopes of the 
surface. This plaquette configuration models a valley bot- 
tom site. 
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Figure 3: Scaling behavior of the isotropic LCOD model for 
L = 16,32,64,128 (top to bottom at the left side). The 
data collapse (/(«)) has been achieved with the MH class 
exponents with conservative noise lll2ll . The insert shows the 
finite size (logarithmic) scaling in the steady state. 



we imposed additional constraints for the accept- 
ing a move. We have tried two kinds of rules based 
on the local neighborhood configurations. The first 
one requires that the height of a particle at the final 
state is higher than that of its initial site 



h&n ~ h ini > 



(15) 



which makes the surface more rough (inverse or 
roughening diffusion). For details of this so-called 
Larger Height Octahedron Model (LHOD) see f^ . 
Note, that LHOD introduces up-down anisotropy, 
which leads to a V 2 [V(/i) 2 ] type of nonlinearity and 
Molecular Beam Epitaxy (MBE) class scaling [29| 
in general. 

The second one is the so-called Larger Curvature 
Octahedron Model (LCOD), which satisfies the de- 
tailed balance condition ([1]), enabling us to realize 
linear, equilibrium MH diffusion steps. The local 
curvature of the surface is calculated at the 4 edges 
of squares of the projected octahedra as shown on 
Fig. [21 We followed the scaling behavior of our 
LCOD in case of smoothing reactions, correspond- 
ing to k > in Eq. 21 Due to the purely diffusive 
noise the growth dynamics is logarithmically slow 
and we find the emergence of the class ([12")) scaling 
(see Figure [3]). 

If the accepted move increases the local curva- 
ture, we call it inverse MH (iMH) process. In the 
continuum model this should correspond to k < 0, 



an unstable growth, however in the finite, lattice 
model formation of pyramid-like structures occurs 
similarly as in the "n = 2" SOS model [30]. When 
we simulated the evolution of the iMH process 
(LCOD model) in the presence of a small, com- 
peting EW type of noise (p = q = 0.05) starting 
from flat initial condition we found a scaling be- 
havior, which agrees well with that of MH univer- 
sality class (fTTj) (see Fig. 0}. The effective growth 
exponent defined as 



/?eff(*) 



In W(t, L -> oo) - In W(t/2, L -> oo) 
ln(f) - ln(i/2) 



(16) 

converges to j3 = 0.25(1) before reaching the satu- 
ration regime. 



2. Pattern generation by competing inverse 
MH and KPZ processes 

A further step in generalizing our basic sur- 
face models was the combination of different sub- 
processes, resulting in noncquilibrium system with 
patterns. For example the mixed application of 
iMH diffusion steps with KPZ updates can model a 
noisy Kuramoto-Sivashinsky (KS) type of equation 



3J,[32j,[33j, the inverse KS (iKS): 



d t h = crV 2 /i + cV 4 /i + A(V/i) 2 



(17) 
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even, blurred and cut into smaller pieces. If we 




Figure 4: Scaling behavior of the isotropic LCOD model 
for L = 32,64,128 (top to bottom at the left side). The 
data collapse (f(u)) has been achieved with the MH class 
exponents 111! . For the growth exponent fitting (dashed 
line) results in ft = 0.25(1). The right insert shows the same 
via the local slopes (IB) . The left insert is a snapshot of 
surface heights of the dot patterns generated by the LHOD 
model with parameters: D± x = D± y = 1 and p = q = 0.1 
at t = 10 4 MCs 



However the signs of couplings are reversed with 
respect to the KS equation \j. This leads to un- 
stable solutions, the long wave behavior of iKS is 
not defined asymptotically in the continuum model. 
However, in the lattice regularized version, for in- 
termediate times we expect dynamical scaling of the 
KS based on symmetry and renormalization argu- 
ments @,[37|. Since the roughening/smoothing sur- 
face moves correspond to phase separation/mixing 
of the lattice gas system it is reasonable to believe 
that both cases are described by the same universal 
fixed point behavior. 

We performed the surface diffusion steps alter- 
nately with the deposition and the removal pro- 
cesses. Thus terms like n x d 4 h and n y d A h with 
y 11 can be modeled. We followed the sur- 



face roughness and calculated PSD as well as the 
wavelength growth of the pattern formation. 

In case of an anisotropic LHOD a competing EW 
process always generates stable ripple patterns, but 
if strong up-down anisotropy (KPZ instead of EW) 
is applied for late times the ripples become un- 
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Figure 5: Snapshot of surface heights of the ripple patterns 
generated by the parameters D x = D y = 0, D- x = D- y = 1 
(anisotropic, iMH) and p = q = 0.05 adsorption/desorption 
at t = 10 4 MCs, in the LHOD model of linear size L = 512. 



1 Where both the V 2 and the V 4 terms have negative 
coefficients and a positive (V/i) 2 can stabilize the solution. 



introduce x/y lattice anisotropy, for example by: 
D x = D y = and D_ x = D^ y = 1 the com- 
peting weak EW or KPZ moves result in ripple 
coarsening (see Fig. [SJ , with a wavelength growth 
A oc i 019 ! 1 ). We calculated this behavior from 
the maxima of the PSD functions. However, for 
late times this coarsening slows down for LCOD 
(A oc t°- 12 W) or becomes faster in case of nonlinear 
LHOD (A oc t°- 35 ( 5 )) as shown on Fig. El The last 
value agrees with that of the numerical estimates 
obtained for the two-field model j38|. 

When isotropic, iMH competes with a weak EW 
process one can observe dot formation in both the 
LHOD and LCOD model. Again, for strong par- 
ticle deposition or removal the patterns fade away 
for long times and the KPZ scaling emerges. The 
insert of Fig. [4] shows a snapshot of the growing 
dots in the LHOD model in the presence of weak 
EW processes. We have determined the PSDs at 
different times for L = 1024 and deduced the time 
dependence of the characteristic wavelength from 
the maxima of 5(k, t) (see Fig. [7]). 

Having confirmed that the LCOD model exhibits 
MH class scaling, we introduced a model for the 
iKS equation (IT?]) a s the combination of iMH and 
KPZ moves. In [29( we presented extensive simu- 
lations providing numerical evidence for an asymp- 
totic KPZ scaling of this model. Note, that unlike 
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Figure 6: PSD for anisotropic diffusion with D y = 1, p = q = 
0.005 for L = 1024 at different times. The insert shows the 
wavelength coarsening calculated from the maxima of the 
PSD curves. Squares: LHOD model with fastening ripple 
formation A oc t0- 35 ( 5 ) ) circles: LCOD model with p = q = 0. 
The lines correspond to power-law fits: A oc iP- la W initially 
and A oc t° 12 ( 1 ) for late times. 



Figure 7: PSD in the LCOD model for isotropic surface 
diffusion (D± x = D± y = 1) with competing EW moves: 
p = q = 0.005 at different times. The insert shows the wave- 
length coarsening calculated from the maxima of the PSD 
curves. Circles correspond to the peaks of the PSD curves, 
dashed line: power- law fit A oc t0.20(l)_ The boxes show the 
same for p = q = in the LCOD. For late times the slope of 
A(i) decreases. 



the KS, the iKS equation exhibits unstable wave- 
length amplification. This can be derived in the 
Fourier space by a linear stability analysis. How- 
ever the dynamical scaling behavior of the lattice 
model corresponding to the iKS equation may be 
expected to be the same as that of the KS equation 
by symmetry arguments. 

3. Conclusions and outlook 

We have returned to some unresolved questions 
of basic surface growth phenomena using mappings 
and extensive computer simulations. Contrary to 
many earlier studies we were able to perform nu- 
merical analysis in 2 + 1 dimensions due to the ef- 
ficient binary lattice gas representation. An im- 
portant advantage of our method as compared to 
the analytical approaches is that we can describe 
anisotropic surface diffusion terms K x d 4 h, K y d 4 h by 
our model. 

We have shown that in the zero external noise 
limit RSOS models can be constructed with short 
range interactions exhibiting MBE or MH type of 
surface growth behavior. For inverse (roughening) 
diffusion, which increases the unstable growth of lo- 
cal curvatures, pyramid-like structures can emerge. 
The size of these structures is limited only by L, 
which is not directly comparable with experiments. 
We created these microscopic models in order to 



study them in competition with the non-conserved 
KPZ processes. 

By mapping surface models onto lattice gases, 
our results imply that the scaling behavior of the 
oriented diffusion of dimers (KPZ) is stable against 
the introduction of an attracting force among them. 
However, a strong repulsion force destroys the fluc- 
tuations, resulting in mean-field class behavior (2^ |. 
We found numerical evidence, by surface scaling 
and probability distribution studies, that the iKS 
model exhibits KPZ scaling in 2 + 1 dimensions. 
The construction of a KS model within the frame- 
work of the lattice gas mapping is under way. Fur- 
ther studies with different boundary conditions can 
reveal interesting connections of the surface tilt to 
the particle concentration of the lattice gas. We 
emphasize that results for disorder dependence or 
anomalous diffusion 36j can easily be transformed 
between these surface and lattice gas models. Ex- 
tension of our mapping can also help to solve more 
detailed models of ion-beam induced pattern forma- 
tion [35 , 34 j , describing the observed narrow band 
of unstable wave-lengths or hexagonal dot forma- 
tion 0. 

We have also determined the time depen- 
dent power spectrum and deduced the wavelength 
growth of the patterns. We have investigated this 
for LHOD and LCOD process, with and without 
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spatial anisotropics. In case of uniaxial surface dif- 
fusion ripple morphologies have been found, while 
for x/y lattice isotropy dot like pattern formation 
could be traced. The wavelength growth is slow for 
LCOD models. Usually we found a X(t) oc t°- w M 
time dependence, which slows down further at late 
times for the linear LCOD model, but becomes 
faster in the nonlinear LHOD model. This agrees 
with the two-field model results [38[. In [29( we 
estimated the characteristic length scale by a dif- 
ferent measure, the average length of the longest 
slopes. The growth of that scale was found to be 
even slower, except when steady DC current flowed 
through the system. Since that quantity measures 
an extremal property of the patterns one can un- 
derstand that somewhat different results arise from 
those of the PSD analysis. 

Finally we point out that these models enable 
efficient, bit-coded, stochastic cellular automaton 
type of simulation of surfaces which can be run ex- 
tremely fast on advanced graphic cards 27, 28| . 
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